/*
 * KineticSolver.h
 *
 *  Created on: 28 Jul 2011
 *      Author: Allan
 */

#ifndef KINETICSOLVER_H_
#define KINETICSOLVER_H_

// C++ includes
#include <string>
#include <vector>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "Core/KineticProblem.h"
#include "Core/Multiphase.h"
#include "Core/ReactionSystem.h"
#include "Core/TotalEquilibriumSolver.h"
#include "Numerics/OdeSolver.h"

class KineticSolver
{
public:
	KineticSolver();
	
	KineticSolver(const Multiphase& multiphase, const ReactionSystem& eReactions, const ReactionSystem& kReactions, 
		const vector<string>& jSpecies,const vector<string>& eSpecies,const vector<string>& kSpecies);
	
	~KineticSolver();
	
	void Initialiaze(const Multiphase& multiphase, const ReactionSystem& eReactions, const ReactionSystem& kReactions, 
		const vector<string>& jSpecies,const vector<string>& eSpecies,const vector<string>& kSpecies);
	
	void SetAlgorithm(OdeAlgorithm algorithm);
	
	void SetControlSolution(double epsAbs, double epsRel);
	
	bool Iterate(double& t, double tEnd, double& dt, double T, double P, VectorXd& n);
	
private:
	/// The multiphase object responsible to compute the activities
	Multiphase multiphase;
	
	/// The equilibrium-controlled reactions
	ReactionSystem eReactions;

	/// The kinetically-controlled reactions
	ReactionSystem kReactions;
	
	/// The kinetics problem, which contains the definition of the ODE function to advance the solution in time 
	KineticProblem problem;
	
	/// The indexes of the primary (jIndexes), equilibrium (eIndexes) and kinetic (kIndexes) species
	vector<Index> jIndexes, eIndexes, kIndexes;
	
	/// The canonical stoichiometric matrix (vej) = -inv(uee) * (uej)
	MatrixXd vej;
	
	/// The number of the primary (Nj), equilibrium (Ne) and kinetic (Nk) species
	unsigned Nj, Ne, Nk;
	
	/// The number of moles of the primary (nj), equilibrium (ne) and kinetic (nk) species
	VectorXd nj, ne, nk;
	
	/// The total number of moles of the primary species (uj) and the vector (ujnk) = [uj nk]
	VectorXd uj, ujnk;
	
	/// The ODE solver
	OdeSolver odeSolver;
	
	/// The total equilibrium solver
	TotalEquilibriumSolver eSolver;
	
	friend class KineticProblem;
};

#endif /* KINETICSOLVER_H_ */
